function eq=SS_eq_CBDCR(func_var,pi_m,i_reserve,i_cbdc,Nb,y)
rq_ratio = func_var.rq_ratio;cost=func_var.cost;
rho_bar= func_var.rho_bar;beta = func_var.beta;df_inv = func_var.df_inv;df = func_var.df;
pi_cbdc = (1+pi_m)./(1+i_cbdc)-1;
pi_reserve_eff = (1+pi_m)./(1+i_reserve)-1;

pi_eff = min(pi_reserve_eff,pi_cbdc);
rho_low = 1/(1+pi_eff)-1;
xi_func=@(x)(1-rq_ratio).*max(1+x,1/(1+pi_eff))+rq_ratio/(1+pi_eff)-cost;
D_bar = max(y(:,1));
psi_func=@(x)interp1(y(:,1),y(:,3),min(x,D_bar)).*(x<=D_bar)+beta.*(x>D_bar);
dpsi_func=@(x)interp1(y(:,1),y(:,4),min(x,D_bar)).*(x<=D_bar);
z_func = @(x)interp1(y(:,1),y(:,2),min(x,D_bar));
Dlow = fzero(@(x)xi_func(rho_low)*(dpsi_func(x).*x/Nb+psi_func(x))-1,[0,10]);
LS_low = (1-rq_ratio)*Dlow*psi_func(Dlow);
LD_high = df_inv(rho_low);
Dbar = fzero(@(x)xi_func(rho_bar).*(dpsi_func(x).*x/Nb+psi_func(x))-1,[0,D_bar - 1e-3]);
LS_bar = (1-rq_ratio)*Dbar*psi_func(Dbar);
LD_low = df_inv(rho_bar);
rho_func = @(x)df(x)-1;
eq_D_func = @(x)xi_func(rho_func(psi_func(x).*x*(1-rq_ratio))).*(dpsi_func(x).*x/Nb+psi_func(x))-1;

if LS_low>LD_high
    eq = [z_func(Dlow),Dlow,rho_low];
elseif LS_bar<LD_low
    eq = [z_func(Dbar),Dbar,rho_bar];
else
    D = fzero(@(x)eq_D_func(x),[Dlow,Dbar]);
    eq = [z_func(D),D,rho_func(D*psi_func(D)*(1-rq_ratio))];
end
    
